AA-GP: process null distribution

This notebook reads in results for the MultiBLUP model using 5000 random gene group to do the following:

• summarize the number of models that pass filtering criteria
  
• examine the distribution of the likelihood ratio (see Edwards et al 2016)
  
• establish 95 percentiles for heritability and the likelihood ratio

knitr::opts_knit$set(root.dir = "~/Documents/AA-GenomicPrediction/",
                     warning = FALSE)
knitr::opts_chunk$set(fig.width = 10)

library(tidyverse)
library(knitr)
library(kableExtra)
library(quantreg)
library(data.table)
library(viridis)
library(emdbook)
library(gridExtra)
library(cowplot)

Size distribution of null pathways

Number of markers should have a uniform distribution

aa_fam <- read.csv("data/processed/aa_family_ids.csv",
                   header = TRUE)

null_dist <- read_csv("reports/lr_null_results.csv") %>%
  # set amino acid family ids 
  left_join(., aa_fam, by = "trait")

null_dist %>%
  distinct(pathway, size.x) %>%
  ggplot(., aes(size.x)) +
  geom_histogram(binwidth = 1000) +
  xlab("Gene group size (number of SNPs)")

Filtering

Filter out values where the REML model did not converge

Set negative heritability estimates to zero

Export summary to S2 table

S2table <- null_dist %>%
  group_by(trait) %>%
  filter(Component %in% "Share_K1") %>%
  mutate("Number of gene groups" = ifelse(SD < 0 | LR < 0, 0, 1),
         "Failed to converge" = ifelse(SD < 0, 1, 0),
         "Negative LR" = ifelse(LR < 0, 1, 0)) %>%
  select("Number of gene groups", "Failed to converge", "Negative LR") %>%
  summarize_all(.funs = "sum") 

Adding missing grouping variables: `trait`

S2table %>%
  write_csv(., "reports/tableS2_null_summary.csv")

head(S2table)

Apply filters to null distribution

null_dist <- null_dist %>% 
  # mutate(Share = ifelse(Share < 0, 0, Share),
  #        Share = ifelse(Share > 1, 1, Share)) %>%
  rename(size = "size.x",
         n_genes = "n_genes.x") %>% 
  # remove values where SD of heritability is strange (results in inflated LR)
  filter(# SD > 0 & LR >= 0 &
    Component %in% "Share_K1") %>%
  dplyr::select(family, trait, pathway, size, n_genes, group_size, LR, Component, Share, SD) 

head(null_dist)

Likelihood ratio

Distribution

Expect LR to be distributed as \(\chi^2_{df 1}\)~ based on Wilks’ Theorem

null_dist %>%
    ggplot(., aes(LR)) + 
    geom_histogram(bins = 50) +
  facet_wrap(~trait, scales = "free", ncol = 6) +
  theme_bw() 

# ggsave("reports/figures/null_histogram.png", height = 12, width = 8)

QQ plots

Maybe some inflation at smaller group sizes (< 10,000 SNPs)

null_dist %>%
  # filter(trait %in% c("ala", "ala_PyrFam", "asp_t", "his", "Total")) %>% 
  mutate(group_size = fct_relevel(group_size, "[0,10000]")) %>%
  ggplot(., aes(sample = LR, colour = group_size)) + 
  scale_colour_manual(values = viridis(6)[1:5]) + 
  stat_qq(distribution = qchisq, dparams = list(df = 1.5)) +
  facet_wrap( ~ trait, scales = "free", ncol = 6) +
  geom_abline(slope = 1, intercept = 0) +
  theme_bw() + 
  theme(legend.position = "top") 

# ggsave("reports/figures/FigS5.png", height = 12, width = 8)

LR goodness-of-fit tests

Test if LR follows Wilks’ theorem expectation of \(\chi^2\) distribution

• \(\chi^2\) test
  
• Kolmogorov-Smirnov test (D)

GoF <- null_dist %>% 
  select(trait, LR) %>%
  group_by(trait) %>%
  mutate(
    # return distance (D) for Kolmogorov-Smirnov goodness-of-fit test 
    # df = 1
    D1 = ks.test(LR, pchisq, df = 1)$statistic,
    # df = mixture of 1 and 2 
    D1_2 = ks.test(LR, rchibarsq(length(LR), 2, mix = 0.5))$statistic,
    # df = 2
    D2 = ks.test(LR, pchisq, df = 2)$statistic,
    # return p-value for chi-sq goodness of fit test 
    # df = 1
    X2_1 = suppressWarnings(chisq.test(LR, rchisq(length(LR), 1))$p.value),
    # df = mixture of 1 and 2
    X2_1_2 = suppressWarnings(chisq.test(LR, rchibarsq(length(LR), 2, mix = 0.5))$p.value),
    # df = 2
    X2_2 = suppressWarnings(chisq.test(LR, rchisq(length(LR), 2))$p.value)) %>%
  summarise_all(mean) %>%
  select(trait, LR, D1, D1_2, D2, X2_1, X2_1_2, X2_2)

# return table of results 
GoF %>% 
  kable(digits = 4) %>%
  kable_styling(bootstrap_options = "striped")

trait	LR	D1	D1_2	D2	X2_1	X2_1_2	X2_2
ala	0.7433	0.0998	0.2528	0.4018	0.2821	0.2821	0.2821
ala_PyrFam	3.6045	0.7322	0.6888	0.6492	0.2741	0.2741	0.2741
ala_t	1.0870	0.0310	0.1540	0.3021	0.2731	0.2731	0.2731
arg	1.4788	0.2762	0.1274	0.1341	0.2631	0.2631	0.2631
arg_GluFam	0.3714	0.2746	0.4298	0.5769	0.2997	0.2997	0.2997
arg_t	0.5028	0.1906	0.3290	0.4820	0.2879	0.2879	0.2879
asp	0.6665	0.1086	0.2638	0.4059	0.2827	0.2827	0.2827
asp_AspFam	1.7296	0.2077	0.0846	0.1587	0.2611	0.2611	0.2611
asp_t	1.7883	0.5250	0.3986	0.2375	0.2666	0.2666	0.2666
AspFam	0.6107	0.1367	0.2860	0.4300	0.2856	0.2856	0.2856
AspFam_Asp	1.6519	0.1948	0.0866	0.1760	0.2620	0.2620	0.2620
BCAA	0.5036	0.1833	0.3398	0.4773	0.2855	0.2855	0.2855
gln	0.4680	0.2112	0.3460	0.5012	0.2872	0.2872	0.2872
gln_GluFam	0.3530	0.3207	0.4626	0.6198	0.3076	0.3076	0.3076
gln_t	0.3872	0.2956	0.4464	0.5948	0.3045	0.3045	0.3045
glu	0.4454	0.2118	0.3526	0.5077	0.2893	0.2893	0.2893
glu_GluFam	0.4530	0.2315	0.3810	0.5324	0.2956	0.2956	0.2956
glu_t	0.9756	0.0514	0.1918	0.3478	0.2776	0.2776	0.2776
GluFam	0.7638	0.1049	0.2222	0.3325	0.2718	0.2718	0.2718
GluFam_glu	0.4788	0.2183	0.3710	0.5199	0.2935	0.2935	0.2935
gly	0.4790	0.2193	0.3564	0.5036	0.2823	0.2823	0.2823
gly_SerFam	0.3489	0.2950	0.4408	0.5970	0.3020	0.3020	0.3020
gly_t	0.2541	0.3801	0.5304	0.6819	0.3093	0.3093	0.3093
his	3.2382	0.5027	0.4780	0.4588	0.2677	0.2677	0.2677
his_GluFam	1.0107	0.0314	0.1876	0.3334	0.2752	0.2752	0.2752
his_t	0.7742	0.0967	0.2296	0.3924	0.2804	0.2804	0.2804
ile	0.5745	0.1624	0.3138	0.4646	0.2866	0.2866	0.2866
ile_AspFam	0.6032	0.1544	0.3024	0.4568	0.2872	0.2872	0.2872
ile_BCAA	3.1959	0.6695	0.6028	0.5313	0.2671	0.2671	0.2671
ile_t	0.8015	0.0837	0.2388	0.3824	0.2814	0.2814	0.2814
leu	0.5778	0.1496	0.2926	0.4426	0.2843	0.2843	0.2843
leu_BCAA	0.7253	0.0989	0.2604	0.3997	0.2830	0.2830	0.2830
leu_PyrFam	1.1848	0.1292	0.1300	0.2752	0.2714	0.2714	0.2714
leu_t	0.9978	0.0251	0.1642	0.3169	0.2748	0.2748	0.2748
lys	0.6401	0.1277	0.2592	0.3825	0.2787	0.2787	0.2787
lys_AspFam	3.0090	0.6554	0.5930	0.5343	0.2718	0.2718	0.2718
lys_t	0.4535	0.2243	0.3788	0.5265	0.2934	0.2934	0.2934
met	0.5487	0.1559	0.2996	0.4417	0.2834	0.2834	0.2834
met_AspFam	-6.3774	0.9996	0.9996	0.9996	0.2861	0.2861	0.2861
met_t	0.5298	0.1675	0.3168	0.4633	0.2865	0.2865	0.2865
phe	0.9022	0.0539	0.2036	0.3560	0.2788	0.2788	0.2788
phe_ShikFam	0.4992	0.1990	0.3444	0.4984	0.2903	0.2903	0.2903
phe_t	1.0603	0.0216	0.1622	0.3115	0.2730	0.2730	0.2730
pro	0.4705	0.1964	0.3382	0.4860	0.2889	0.2889	0.2889
pro_GluFam	0.5255	0.1719	0.2962	0.4496	0.2829	0.2829	0.2829
pro_t	0.4862	0.1919	0.3382	0.4794	0.2871	0.2871	0.2871
PyrFam	0.5045	0.1926	0.3556	0.4928	0.2909	0.2909	0.2909
ser	0.4128	0.2313	0.3884	0.5297	0.2930	0.2930	0.2930
ser_SerFam	0.4666	0.2395	0.3872	0.5395	0.2989	0.2989	0.2989
ser_t	0.4228	0.2611	0.4080	0.5560	0.3015	0.3015	0.3015
SerFam	0.5007	0.2011	0.3414	0.4788	0.2827	0.2827	0.2827
ShikFam	1.6064	0.1484	0.0436	0.1879	0.2628	0.2628	0.2628
thr	1.4055	0.1325	0.0666	0.2166	0.2657	0.2657	0.2657
thr_AspFam	0.7956	0.0949	0.2530	0.3914	0.2817	0.2817	0.2817
total	0.5574	0.1456	0.2850	0.4249	0.2818	0.2818	0.2818
trp	0.6488	0.1471	0.2976	0.4487	0.2861	0.2861	0.2861
trp_ShikFam	0.4828	0.2035	0.3354	0.4928	0.2859	0.2859	0.2859
trp_t	0.5313	0.1832	0.3136	0.4585	0.2813	0.2813	0.2813
tyr	0.5001	0.1823	0.3240	0.4780	0.2882	0.2882	0.2882
tyr_ShikFam	0.4890	0.1888	0.3278	0.4752	0.2862	0.2862	0.2862
tyr_t	0.8177	0.0762	0.2144	0.3606	0.2811	0.2811	0.2811
val	0.6674	0.1248	0.2408	0.3726	0.2748	0.2748	0.2748
val_BCAA	1.3624	0.4419	0.2864	0.2004	0.2746	0.2746	0.2746
val_PyrFam	0.6507	0.1195	0.2688	0.4137	0.2829	0.2829	0.2829
val_t	0.3689	0.2647	0.4112	0.5667	0.2959	0.2959	0.2959

# plot results
GoF %>%
  select(trait, D1, D1_2, D2) %>%
  gather(key = df, value = D, -trait) %>%
  mutate(df = recode(df, D1 = "1", D1_2 = "1.5", D2 = "2")) %>%
  ggplot(., aes(df, D, group = trait)) +
  geom_point() + 
  geom_line() +
  xlab("degrees of freedom") +
  ylab("Kolmogorov-Smirnov test statistic (D)") + 
  facet_wrap(~ trait, scales = "fixed", ncol = 6) +
  theme_bw() +
  theme(legend.position = "none")

# ggsave("reports/figures/FigS4.png", height = 12, width = 10)  

GoF %>%
  select(trait, D1, D1_2, D2) %>%
  gather(key = df, value = D, -trait) %>%
  group_by(trait) %>%
  mutate(df = recode(df, D1 = "1", D1_2 = "1.5", D2 = "2")) %>%
  slice(which.min(D)) %>%
  ungroup() %>%
  write_csv("reports/goodness_of_fit.csv") %>%
  head()

Empirical thresholds

95 percentile for LR

Fit an additive quantile regression to establish 95 percentile cut off for LR at different pathway sizes (number of markers)

## 95% LR threshold for plotting
lr_fit <- NULL
coefs <- NULL

for (i in unique(null_dist$trait)) {
  lr_fit[[i]] <- rqss(as.numeric(LR) ~ qss(size, constraint = "I"), 
                      data = null_dist[null_dist$trait==i,], tau = 0.95) 
  coefs[[i]] <- data.frame(lr_fit[[i]]$qss$size$xyz) %>%
    mutate(lr_95 = X2 + lr_fit[[i]]$coef[1],
           size = X1) %>%
    dplyr::select(size, lr_95)
}

lr_threshold <- rbindlist(coefs, idcol = "trait")

null_dist <- null_dist %>%
  left_join(., lr_threshold, by = c("trait", "size")) %>%
  mutate(col = LR >= lr_95) 

null_dist %>%
  as_tibble() %>%
  head()

95 percentile for \(h^2\)

Fit an additive quantile regression to establish 95% quantile cut off for proportion of genomic variance explained (H^2) at different pathway sizes (number of markers)

h2_fit <- NULL
h2_coefs <- NULL

for (i in unique(null_dist$trait)) {
  h2_fit[[i]] <- rqss(as.numeric(Share) ~ qss(size, constraint = "I"), 
                      data = null_dist[null_dist$trait==i,], tau = 0.977) #FDR corrected 
                      # data = null_dist[null_dist$trait==i,], tau = 0.95)
  h2_coefs[[i]] <- data.frame(h2_fit[[i]]$qss$size$xyz) %>%
    mutate(h2_95 = X2 + h2_fit[[i]]$coef[1],
           size = X1) %>%
    dplyr::select(size, h2_95)
}

h2_threshold <- rbindlist(h2_coefs, idcol = "trait") %>%
  as_tibble()
  
head(h2_threshold)

Plot null distribution

Blue dots are pathways that also pass the LR_95 threshold
Solid line is the 95% percentile for proportion of heritability explained
Dashed line is the infinitesimal expectation (each SNP contributes approximately equal amount of variation)

save results

save(lr_fit, h2_fit, file = "reports/null_distribution.RData")

Error in save(lr_fit, h2_fit, file = "reports/null_distribution.RData") : 
  object ‘h2_fit’ not found